sets             j       functions                       /1*3/
                 e       finite elements                 /1*4/
                 c       collocation points              /1*4/
                 i       unknown parameters              /1*3/
                 m       measurements                    /1*3/
                 p       patients                        /1*40/
                 em      element associated with measurement;
alias(c,cc,ccc,c4);
alias(m,mm);


parameters       initial_condition
                 f_measured      value of measurement m for function j
                 t               time of measurement m
                 tau             normalized legendre collocation points
                 tau_m           normalized time points for measurements
                 ti              initial time of element e
                 tf              final time of element e;

Scalar
         pi                      pi number               /3.14159265/
         precision_posterior                             /0.25/;


Table data_Y2(*,p)
$ondelim
$include data_GAMS.csv
$offdelim
;

Table time(*,p)
$ondelim
$include time_GAMS.csv
$offdelim
;
Table data_patients(p,p)
$ondelim
$include initial_conditions.csv
$offdelim
;

variables        param           parameters to be estimated
                 f_value         function value at the colloc pt c in elem e for func j
                 error
                 log_likelihood;

positive variable   f_predicted     predicted function value at observed points;

equations        residuals1      residual for the differential-algebraic part
                 residuals2
                 continuity      continuity profile equations for the elements
                 predictions     equation to compute predictions for measurement
                 mass_balance    mass balance
                 error_eq
                 OF              objective function;

initial_condition('1',p)=data_patients(p,'1');
initial_condition('2',p)=data_patients(p,'2');

tau('1')=-1;
tau('2')=-0.774596669241;
tau('3')=0;
tau('4')=0.774596669241;


f_measured('1',m,p)=0;
f_measured('2',m,p)=data_Y2(m,p);
t(m,p)=time(m,p);


ti('1')=0;
tf('1')=27;
ti('2')=27;
tf('2')=31;
ti('3')=31;
tf('3')=109;
ti('4')=109;
tf('4')=180;

em(e,m,p)=no;
em(e,m,p)$(t(m,p) ge ti(e) and t(m,p) lt tf(e))=yes;
em(e,m,p)$(t(m,p) eq tf(e) and ord(e) eq card(e))=yes;

loop((e,m,p)$em(e,m,p),
         tau_m(m,p)=1-2*(tf(e)-t(m,p))/(tf(e)-ti(e));
);


display t,ti,tf,em,tau_m;


residuals1(j,e,c,p)$(ord(c) gt 1 and ord(j) eq 1)..
         sum(cc,
                 f_value(j,e,cc,p)*sum(c4$(ord(c4) ne ord(cc)), prod(ccc$(ord(ccc) ne ord(cc) and ord(ccc) ne ord(c4)), (tau(c)-tau(ccc))))/
                 prod(ccc$(ord(ccc) ne ord(cc)), (tau(cc)-tau(ccc)))
         )=e=(tf(e)-ti(e))*((-1)*param('1')*f_value('1',e,c,p)*f_value('3',e,c,p))/2;

residuals2(j,e,c,p)$(ord(c) gt 1 and ord(j) eq 2)..
         sum(cc,
                 f_value(j,e,cc,p)*sum(c4$(ord(c4) ne ord(cc)), prod(ccc$(ord(ccc) ne ord(cc) and ord(ccc) ne ord(c4)), (tau(c)-tau(ccc))))/
                 prod(ccc$(ord(ccc) ne ord(cc)), (tau(cc)-tau(ccc)))
         )=e=(tf(e)-ti(e))*(param('2')*f_value('3',e,c,p))/2;

continuity(j,e,p)$(ord(j) le 2)..
         f_value(j,e,'1',p)=e=sum(c,
                 f_value(j,e-1,c,p)*
                 prod(cc$(ord(cc) ne ord(c)), (1-tau(cc))/(tau(c)-tau(cc)))
         )+initial_condition(j,p)$(ord(e) eq 1);


predictions(j,m,p)$(ord(j) eq 2)..
         f_predicted(j,m,p)=e=sum((e,c)$em(e,m,p),
                 f_value(j,e,c,p)*
                 prod(cc$(ord(cc) ne ord(c)), (tau_m(m,p)-tau(cc))/(tau(c)-tau(cc)))
         );

mass_balance(e,c,p)..  f_value('1',e,c,p)+f_value('2',e,c,p)+f_value('3',e,c,p)=e=100;

error_eq(j,m,p)$(ord(j) eq 2)..
              error(j,m,p)=e=log(f_measured(j,m,p))-log(f_predicted(j,m,p));

OF..
         log_likelihood=e=-sum((j,m,p)$(ord(j) eq 2),
                 0.5*log(param('3')/(2*pi))-0.5*param('3')*error(j,m,p)*error(j,m,p)
                 );

param.lo('1')=0;
param.lo('2')=0;
param.up('1')=25;
param.up('2')=25;
param.lo('3')=0.000001;
param.up('3')=100000;
f_predicted.lo(j,m,p)=0.000000000001;

param.fx('1')=0.0005;
param.fx('2')=8;
param.fx('3')=1;

option limrow=1e3;
option optca=1.5e-2;
option optcr=0.03;
option decimals=6;
model Var_Infer /all/;


option nlp=conopt;
Var_Infer.optfile=1
solve Var_Infer using NLP min log_likelihood;
Execute_Unload 'PhaseI';

display param.l, log_likelihood.l;


*$ontext
option nlp=conopt;
param.lo('1')=0;
param.lo('2')=0;
param.up('1')=25;
param.up('2')=25;
param.lo('3')=0.000001;
param.up('3')=100000;
f_predicted.lo(j,m,p)=0.000000000001;
Execute_loadpoint 'PhaseI';
solve Var_Infer using NLP min log_likelihood;
Execute_Unload 'PhaseI';
*$offtext

$ontext
option nlp=conopt;
param.lo('1')=0.00001;
param.lo('2')=0.00001;
param.up('1')=25;
param.up('2')=25;
param.lo('3')=0.000001;
param.up('3')=10000;
*f_predicted.lo('2',m,p)=0.5*f_measured('2',m,p);
f_predicted.up('2',m,p)=1.5*f_measured('2',m,p);
*f_value.lo('2',e,c,p)=0.0*smin((m), f_measured('2',m,p));
f_value.up('2',e,c,p)=1.5*smax((m), f_measured('2',m,p));
Execute_loadpoint 'PhaseI';
solve Var_Infer using NLP min log_likelihood;
$offtext




parameter CPUsec1, meanPhase1, max_posterior;
FILE phase1 /phase1.csv/;
phase1.PC=5;
phase1.nd=7;
phase1.pw=2000;
PUT phase1;
CPUsec1=Var_Infer.resusd;
meanPhase1(i)=log(param.l(i));
meanPhase1('3')=log(sqrt(1/param.l('3')));
max_posterior=-log_likelihood.l;
put CPUsec1;
loop(i,
         put meanPhase1(i);
);
put max_posterior;
putclose;

display param.l, f_measured, f_predicted.l, f_value.l, log_likelihood.l, meanPhase1, error.l;



